##################################################################
# File:      contrasts_online_appendix.R
# Purporse:  Replication of the online appendix figures and tables for 
#                 Journal of Conflict Resolution
# Author:    Nadiya Kostyuk
# R version 4.2.2 (2022-10-31 ucrt)
# Last Edit: 01/05/2023
##################################################################

# model fitting
# contrasts and regression table
# figures and plots

rm(list=ls())

## Install & load packages (all at once)
list.of.packages <- c('dplyr', 'doBy')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}
lapply(list.of.packages, require, character.only = TRUE)

# setting directory: 
setwd("C:/Users/nkostyuk3/Dropbox (GaTech)/ResearchProjects/GartzkeKostyuk/Complements v. Substitutes/Publication/Complementarity/JCR/Replication")

#########################################################
# Example of how to estimate contrasts without covariates
#########################################################

# loading the model results obtained earlier:
load('mod_list_tab2a.RData')

first_mod = seq(1:6)
i = 1; first_mod[i]

results_list1 = list()

for(i in 1:length(first_mod)){
  
  ######
  # first model (no covariates)
  L1 <- matrix(0, nrow=6, ncol = length(summary(mod_list[[i]])$coefficients[,1]))
  L1
  colnames(L1) = names(summary(mod_list[[i]])$coefficients[,1])
  rownames(L1) = c('P_cyber_base', 'P_military_base',
                   'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                   'OR_military_mil(lag)', 'OR_cyber_mil(lag)') # OR (odds ratio)
  L1[1,c('(Intercept)', 'TYPE')] <- c(1,1) 
  L1[2,c('(Intercept)', 'TYPE')] <- c(1,-1) 
  L1[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1) 
  L1[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
  L1[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
  L1[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
  L1
  res1 <- doBy::esticon(mod_list[[i]], L=L1)
  
  res1 = res1 %>%
    dplyr::mutate(starz=case_when(
      p.value>=.1~'',
      p.value>=0.05&p.value<0.1~'^',
      p.value>=0.01&p.value<0.05~'*',
      p.value>=0.001&p.value<0.01~'**',
      p.value<0.001~'***'
    ))
  res1
  
  m=1.96
  res1_1 <- with(res1,
                 sprintf('%4.2f%s(%4.2f;%4.2f)',
                         exp(estimate),
                         starz,
                         exp(estimate-m*std.error),
                         exp(estimate+m*std.error)
                 ) 
  ) 
  res1_1
  names(res1_1) = rownames(L1)
  res1_1
  
  aic = round(AIC(mod_list[[i]]), 1)
  nobz = nobs(mod_list[[i]])
  model_name = paste0(srcz[u], lvz[l])
  res1_1 = c(res1_1, aic, nobz, model_name)

  results_list1[[i]] = res1_1
  cat(i, 'out of', length(first_mod), '\n')
}
results_list1
data_rez = data.frame(matrix(unlist(results_list1), nrow=length(results_list1), byrow=TRUE))
colnames(data_rez) = names(results_list1[[1]])
data_rez = t(data_rez)
data_rez
xtable(data_rez)

writeLines(capture.output(xtable(data_rez)), 'table_no_control.tex')


#########################################################
# Example of how to estimate contrasts without covariates
#########################################################

# Note: this is the main model with interation effects for Attacker's (A) and Target's (B) internet access 
# and Attacker's CINC score


L <- matrix(0, nrow=17, ncol = length(summary(mod3)$coefficients[,1]))
L
colnames(L) = names(summary(mod3)$coefficients[,1])
rownames(L) = c('P_cyber_base', 'P_military_base',
                'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                'OR_military_mil(lag)', 'OR_cyber_mil(lag)',
                'OR_conflict_intus(a)','OR_conflict_intus(b)',
                'OR_conflict_cinc(a)', 'OR_conflict_distance',
                'OR_conflict_nuclear(b)',
                'OR_cyber_intus(a)','OR_military_intus(a)',
                'OR_cyber_intus(b)','OR_military_intus(b)',
                'OR_cyber_cinc(a)', 'OR_military_cinc(a)') 
L[1,c('(Intercept)', 'TYPE')] <- c(1,1) 
L[2,c('(Intercept)', 'TYPE')] <- c(1,-1) 
L[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1) 
L[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
L[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
L[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
L[7,c('LOG_INT_USERS_A_sc')] <- 1
L[8,c('LOG_INT_USERS_B_sc')] <- 1
L[9,c('CINC_A_sc')] <- 1
L[10,c('INV_CAP_DISTANCE_sc')] <- 1 
L[11,c('NUCLEAR')] <- 1
L[12,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,1) 
L[13,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,-1)
L[14,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,1) 
L[15,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,-1)
L[16,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,1) 
L[17,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,-1)
L
res <- doBy::esticon(mod3, L=L)
res

res = res %>%
  dplyr::mutate(starz=case_when(
    p.value>=.1~'',
    p.value>=0.05&p.value<0.1~'^',
    p.value>=0.01&p.value<0.05~'*',
    p.value>=0.001&p.value<0.01~'**',
    p.value<0.001~'***'
  ))
res

m=1.96
res3_1 <- with(res[1:2,],
               sprintf('%4.2f (%4.2f, %4.2f)',
                       plogis(estimate),
                       plogis(estimate-m*std.error),
                       plogis(estimate+m*std.error)
               ) # replace prob with 2 decimal places
               
) 
res3_1
names(res3_1) = rownames(L[1:2,])
res3_1

res3_2 <- with(res[3:17,],
               sprintf('%4.2f%s(%4.2f;%4.2f)',
                       exp(estimate),
                       starz,
                       exp(estimate-m*std.error),
                       exp(estimate+m*std.error)
               ) 
               
) 
res3_2
names(res3_2) = rownames(L[3:17,])
res3_2


aic = round(AIC(mod3), 1)
nobz = nobs(mod3)
res3_2 = c(res3_2, aic = aic, nobz = nobz)

writeLines(capture.output(xtable(res3_2)), 'table_main_model.tex')


